#!/usr/bin/env R
# -*- coding: utf-8  -*-

# This file is part of PCR efficiency calculator.  
# PCR efficiency calculator is free software: you can
# redistribute it and/or modify it under the terms of the GNU General Public
# License as published by the Free Software Foundation, version 2.
#
# This program is distributed in the hope that it will be useful, but WITHOUT
# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
# FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more
# details.
#
# You should have received a copy of the GNU General Public License along with
# this program; if not, write to the Free Software Foundation, Inc., 51
# Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
#
# Copyright Izaskun Mallona
# izaskun.mallona@gmail.com

library(Cairo)
library(ROCR)
library(verification)

#cleaning up the predictions
predictions<-round(predict(miniGaussian),3)
predictions[predictions<1.5]<-1.5
predictions[predictions>=2]<-2

#starting verification

ROC<-list(predictions=as.vector(predictions),observations=dataMini$efficiency[!is.na(dataMini$efficiency)])
ROC.60<-list(predictions=as.vector(predictions),observations=dataMini$efficiency[!is.na(dataMini$efficiency)])
ROC.70<-list(predictions=as.vector(predictions),observations=dataMini$efficiency[!is.na(dataMini$efficiency)])
ROC.80<-list(predictions=as.vector(predictions),observations=dataMini$efficiency[!is.na(dataMini$efficiency)])
ROC.90<-list(predictions=as.vector(predictions),observations=dataMini$efficiency[!is.na(dataMini$efficiency)])

#observations will be cut with a threshold of 1.60 of efficiency

ROC.60$observations[ROC$observations<1.60]<-0
ROC.60$observations[ROC$observations>=1.60]<-1

ROC.70$observations[ROC$observations<1.70]<-0
ROC.70$observations[ROC$observations>=1.70]<-1

ROC.80$observations[ROC$observations<1.80]<-0
ROC.80$observations[ROC$observations>=1.80]<-1

ROC.90$observations[ROC$observations<1.90]<-0
ROC.90$observations[ROC$observations>=1.90]<-1

ROCdf<-data.frame()
CairoSVG('verification.svg')

par(mfrow=c(3,3))
## computing a simple ROC curve (x-axis: fpr, y-axis: tpr)
pred <- prediction( ROC.60$predictions, ROC.60$observations)
perf <- performance(pred,"tpr","fpr")
plot(perf,col='blue',type='l',pch=1,xlim=c(0,1),ylim=c(0,1))
## precision/recall curve (x-axis: recall, y-axis: precision)
perf1 <- performance(pred, "prec", "rec")
plot(perf1,col='blue',type='l',pch=1,xlim=c(0,1),ylim=c(0,1))
## sensitivity/specificity curve (x-axis: specificity,
## y-axis: sensitivity)
#perf1 <- performance(pred, "sens", "spec")
#plot(perf1)


## computing a simple ROC curve (x-axis: fpr, y-axis: tpr)
pred <- prediction( ROC.70$predictions, ROC.70$observations)
perf <- performance(pred,"tpr","fpr")
plot(perf,col='red',type='l',pch=2,xlim=c(0,1),ylim=c(0,1))
## precision/recall curve (x-axis: recall, y-axis: precision)
perf1 <- performance(pred, "prec", "rec")
plot(perf1,col='red',type='l',pch=2)
## sensitivity/specificity curve (x-axis: specificity,
## y-axis: sensitivity)
#perf1 <- performance(pred, "sens", "spec")
#plot(perf1)


## computing a simple ROC curve (x-axis: fpr, y-axis: tpr)
pred <- prediction( ROC.80$predictions, ROC.80$observations)
perf <- performance(pred,"tpr","fpr")
plot(perf,col='green',type='l',pch=3,xlim=c(0,1),ylim=c(0,1))
## precision/recall curve (x-axis: recall, y-axis: precision)
perf1 <- performance(pred, "prec", "rec")
plot(perf1,col='green',type='l',pch=3,xlim=c(0,1),ylim=c(0,1))
## sensitivity/specificity curve (x-axis: specificity,
## y-axis: sensitivity)
#perf1 <- performance(pred, "sens", "spec")
#plot(perf1)


## computing a simple ROC curve (x-axis: fpr, y-axis: tpr)
pred <- prediction( ROC.90$predictions, ROC.90$observations)
perf <- performance(pred,"tpr","fpr")
plot(perf,col='black',type='l',pch=4,xlim=c(0,1),ylim=c(0,1))
## precision/recall curve (x-axis: recall, y-axis: precision)
perf1 <- performance(pred, "prec", "rec")
plot(perf1,col='black',type='l',pch=4,xlim=c(0,1),ylim=c(0,1))
## sensitivity/specificity curve (x-axis: specificity,
## y-axis: sensitivity)
#perf1 <- performance(pred, "sens", "spec")
#plot(perf1)

dev.off()

